;;----------------------------------------------------------------------
;; plot zonal distribution of ion production 
;;----------------------------------------------------------------------
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$KAILIBROOT/lib_radon.ncl"

begin

  plotname = "fig_surf" 
  plotform = "eps" 

  np = 9 

  str = (/"","","","","","","","",""/) 

  tiYAxisString = str
  tiYAxisString(0) = "Height above surface (m)" 
  tiYAxisString(3) = "Height above surface (m)" 
  tiYAxisString(6) = "Height above surface (m)" 

  abcd = (/"a","b","c","d","e","f","g","h","i","j","k","l","m","n"/) 

  abcd = abcd + ") " 

  gsnLeftString = str

  stra = "WCRP1995 " 
  strb = "SW1998 / 1.6 " 
  strc = "MERGE " 
 
  gsnLeftString(0) = abcd(0) + stra
  gsnLeftString(1) = abcd(1) + strb
  gsnLeftString(2) = abcd(2) + strc

  gsnLeftString(3) = abcd(3) + stra 
  gsnLeftString(4) = abcd(4) + strb
  gsnLeftString(5) = abcd(5) + strc

  gsnLeftString(6) = abcd(6) + stra 
  gsnLeftString(7) = abcd(7) + strb
  gsnLeftString(8) = abcd(8) + strc 
 
  gsnCenterString = str

  strx = "IPR ( cm~S~-3~N~s~S~-1~N~ ) "
  stra = strx ;;+ "WCRP1995 " 
  strb = strx ;;+ "SW1998 / 1.6 " 
  strc = strx ;;+ "MERGE " 

  gsnCenterString(0) = stra
  gsnCenterString(1) = strb
  gsnCenterString(2) = strc 

  gsnCenterString(3) = stra 
  gsnCenterString(4) = strb
  gsnCenterString(5) = strc 

  gsnCenterString(6) = stra 
  gsnCenterString(7) = strb
  gsnCenterString(8) = strc 

;;----------------------------------------------------------------------
;; load files 
;;----------------------------------------------------------------------

  overland = True
 
  path = "/pf/m/m222044/libdata/"
  pata = "/pf/m/m222044/rnpaper/"


  run1 = "RA008" 
  run2 = "RA010" 
  run3 = "RA007" 

  filenama = "out_"+run1+"_ANN.nc" 
  filenamb = "out_"+run1+"_DJF.nc" 
  filenamc = "out_"+run1+"_JJA.nc" 

  filenamd = path+"mbase_T63L31.nc" 

  filenamf = "out_"+run2+"_ANN.nc" 
  filenamg = "out_"+run2+"_DJF.nc" 
  filenamh = "out_"+run2+"_JJA.nc" 

  filenami = "out_"+run3+"_ANN.nc" 
  filenamj = "out_"+run3+"_DJF.nc" 
  filenamk = "out_"+run3+"_JJA.nc" 


  fila = addfile(filenama,"r")
  filb = addfile(filenamb,"r")
  filc = addfile(filenamc,"r")
  fild = addfile(filenamd,"r")

  filf = addfile(filenamf,"r")
  filg = addfile(filenamg,"r")
  filh = addfile(filenamh,"r")

  fili = addfile(filenami,"r")
  filj = addfile(filenamj,"r")
  filk = addfile(filenamk,"r")

;;----------------------------------------------------------------------
;; read data 
;;----------------------------------------------------------------------

  rna = fila->RN222    ;; radon  
  rnb = filb->RN222
  rnc = filc->RN222

  rma = fila->RN22290  ;; 90th percentile   
  rmb = filb->RN22290
  rmc = filc->RN22290

  roa = fila->RN22250  ;; 90th percentile   
  rob = filb->RN22250
  roc = filc->RN22250 

  rnf = filf->RN222
  rmf = filf->RN22290 

  rng = filg->RN222
  rmg = filg->RN22290 

  rnh = filh->RN222
  rmh = filh->RN22290 

  rni = fili->RN222
  rmi = fili->RN22290 

  rnj = filj->RN222
  rmj = filj->RN22290 

  rnk = filk->RN222
  rmk = filk->RN22290 

  slm = fild->slm

;;----------------------------------------------------------------------
;; new data structure for zonal average 
;;----------------------------------------------------------------------

  nlev = dimsizes(rna(0,:,0,0)) 

  plev = nlev-1
 
  prna = rna(0,plev,:,:) 
  prnb = rnb(0,plev,:,:) 
  prnc = rnc(0,plev,:,:)

  prnf = rnf(0,plev,:,:) 
  prng = rng(0,plev,:,:) 
  prnh = rnh(0,plev,:,:) 

  prni = rni(0,plev,:,:) 
  prnj = rnj(0,plev,:,:) 
  prnk = rnk(0,plev,:,:) 

  prma = rma(0,plev,:,:)
  prmb = rmb(0,plev,:,:)
  prmc = rmc(0,plev,:,:) 

  prmf = rmf(0,plev,:,:) 
  prmg = rmg(0,plev,:,:) 
  prmh = rmh(0,plev,:,:) 

  prmi = rmi(0,plev,:,:) 
  prmj = rmj(0,plev,:,:) 
  prmk = rmk(0,plev,:,:) 

  proa = roa(0,plev,:,:)
  prob = rob(0,plev,:,:)
  proc = roc(0,plev,:,:) 

  prna = RN222_ionprod(prna) 
  prnb = RN222_ionprod(prnb) 
  prnc = RN222_ionprod(prnc) 
 
  prnf = RN222_ionprod(prnf) 
  prng = RN222_ionprod(prng) 
  prnh = RN222_ionprod(prnh) 

  prni = RN222_ionprod(prni) 
  prnj = RN222_ionprod(prnj) 
  prnk = RN222_ionprod(prnk) 

  prma = RN222_ionprod(prma) 
  prmb = RN222_ionprod(prmb) 
  prmc = RN222_ionprod(prmc) 

  proa = RN222_ionprod(proa) 
  prob = RN222_ionprod(prob) 
  proc = RN222_ionprod(proc) 

  prmf = RN222_ionprod(prmf) 
  prmg = RN222_ionprod(prmg) 
  prmh = RN222_ionprod(prmh) 

  prmi = RN222_ionprod(prmi) 
  prmj = RN222_ionprod(prmj) 
  prmk = RN222_ionprod(prmk) 

  if (overland) then 

  pslm = slm(0,:,:) 

  prna = where(pslm.lt.0.5,-999.,prna)  
  prnb = where(pslm.lt.0.5,-999.,prnb)  
  prnc = where(pslm.lt.0.5,-999.,prnc)  

  prnf = where(pslm.lt.0.5,-999.,prnf)  
  prng = where(pslm.lt.0.5,-999.,prng)  
  prnh = where(pslm.lt.0.5,-999.,prnh)  

  prni = where(pslm.lt.0.5,-999.,prni)  
  prnj = where(pslm.lt.0.5,-999.,prnj)  
  prnk = where(pslm.lt.0.5,-999.,prnk)  

  prma = where(pslm.lt.0.5,-999.,prma)  
  prmb = where(pslm.lt.0.5,-999.,prma)  
  prmc = where(pslm.lt.0.5,-999.,prma)   

  prmf = where(pslm.lt.0.5,-999.,prmf)  
  prmg = where(pslm.lt.0.5,-999.,prmg)  
  prmh = where(pslm.lt.0.5,-999.,prmh)  

  prmi = where(pslm.lt.0.5,-999.,prmi)  
  prmj = where(pslm.lt.0.5,-999.,prmj)  
  prmk = where(pslm.lt.0.5,-999.,prmk)  

  end if 

;;----------------------------------------------------------------------
;; open wks  
;;----------------------------------------------------------------------

  wks = gsn_open_wks(plotform,plotname)               ; open a ps file

;;----------------------------------------------------------------------
;; colormap
;;----------------------------------------------------------------------

  colormap = "amwg"
  colormap = "ViBlGrWhYeOrRe"
  colormap = "testcmap"

  gsn_define_colormap(wks,colormap)         ; choose a colormap

;;----------------------------------------------------------------------
;; plot settings:  
;;----------------------------------------------------------------------

  ;; for all panels 
 
  res                  = True    ; contour mods desired
  res@gsnFrame         = False   ; don't draw yet
  res@gsnDraw          = False   ; don't advance frame yet
  res@tiMainOn         = False   ; no title
 ;res@tiYAxisString    = "Height (m)" ; y-axis title
  res@cnFillOn         = True    ; turn on color
  res@cnLinesOn        = False   ; turn off contour lines
  res@cnLineThicknessF = 1.     ; thicker lines
  res@cnLevelSelectionMode  = "ManualLevels"  ; set manual cn levels
  res@cnMinLevelValF   =    0.     ; set min contour level
  res@cnMaxLevelValF   = 8000.     ; set max contour level
  res@cnLevelSpacingF  = 1000.
  res@cnInfoLabelOn    = False    ; no contour labels
 ;res@cnFillMode       = "RasterFill"  ; turn on raster mode
 ;res@cnFillDrawOrder  = "PreDraw"     ; draw contours before continents

  res@gsnSpreadColors = True    ; use full colormap

  if (colormap.eq."testcmap") then 
     res@gsnSpreadColorStart = 50           ; start at color 2
     res@gsnSpreadColorEnd   = 180          ; don't use added gray
  end if 
  if (colormap.eq."ViBlGrWhYeOrRe") then
     res@gsnSpreadColorStart = 24           ; start at color 2
     res@gsnSpreadColorEnd   = 95          ; don't use added gray
  end if
 
  res@gsnLeftString  = "IPR ( cm~S~-3~N~s~S~-1~N~ )"       ;  -~S~222~N~Rn
  res@gsnRightString = ""

 
;;----------------------------------------------------------------------
;; lable bar 
;;----------------------------------------------------------------------
 
  res@lbLabelBarOn       = True       ; no individual label bar
  res@lbLabelStride      = 2          ; every other label bar label
  res@lbOrientation      = "Vertical" ; vertical label bar
  res@lbLabelFontHeightF = 0.018      ; default is HUGE
  res@pmLabelBarWidthF   = 0.08        ; default is shorter
  res@pmLabelBarHeightF  = 0.5        ; default is taller

  res@trYReverse         = False      ; reverse y axis
 ;res@trYMinF            = 0.         ; set minimum Y-axis value
 ;res@trYMaxF            = 8000.      ; set maximum Y-axis value
 ;res@trXMinF            = 1949.      ; set minimum X-axis value
 ;res@trXMaxF            = 2006.      ; set maximum X-axis value

  res@vpWidthF           = 1.0        ; width of contour plots

;;----------------------------------------------------------------------
;; plot array 
;;----------------------------------------------------------------------

  plot = new(6,graphic)                          ; create graphical array


;;----------------------------------------------------------------------
;; ANN panel 
;;----------------------------------------------------------------------

  res@gsnRightString = "ANN"
  res@cnLinesOn       = False    ; turn off contour lines
  res@cnFillOn        = True     ; turn on color
  res@cnMinLevelValF  =    0.    ; set min contour level
  res@cnMaxLevelValF  =   16.      ; set max contour level
  res@cnLevelSpacingF = 2.
  res@lbLabelStride   = 2        ; every other label bar label
  res@gsnLeftString  = "IPR ( cm~S~-3~N~s~S~-1~N~ ) mean"       ;  -~S~222~N~Rn
  plot(0) = gsn_csm_contour_map_ce(wks,prna,res) 

;;----------------------------------------------------------------------
;; line over color 
;;----------------------------------------------------------------------

  res@cnLinesOn       = False    ; turn off contour lines
  res@cnFillOn        = True     ; turn on color
  res@cnLineThicknessF= 2.      ; thicker lines
  res@cnLineLabelFontHeightF = 0.016
 ;res@cnMonoLineDashPattern  = True
  res@cnLineLabelPlacementMode = "Constant" 
  res@cnLineDashPatterns  = 2   ; Define line patterns.
  res@cnLineDashSegLenF   = 0.18          ; assist in controlling concentration
  res@cnLineLabelInterval = 1             ; default = 2 
  res@cnMinLevelValF  =    0.      ; set min contour level
  res@cnMaxLevelValF  =   16.      ; set max contour level
  res@cnLevelSpacingF =    2.
  res@lbLabelStride   =    2       ; every other label bar label
  res@gsnLeftString  = "IPR ( cm~S~-3~N~s~S~-1~N~ ) 90~S~th~N~ percentile"       ;  -~S~222~N~Rn
  plot(1) = gsn_csm_contour_map_ce(wks,prma,res)           ; draw second plot
 ;overlay(plot(0),plotover)                              ; now over lay plots 

 ;delete(plotover) 

;;----------------------------------------------------------------------
;; 2nd panel 
;;----------------------------------------------------------------------

  res@gsnRightString  = "DJF"
  res@cnLinesOn       = False    ; turn off contour lines
  res@cnFillOn        = True     ; turn on color
  res@cnMinLevelValF  =    0.    ; set min contour level
  res@cnMaxLevelValF  =   16.      ; set max contour level
  res@cnLevelSpacingF = 2.
  res@lbLabelStride   = 2        ; every other label bar label
 ;res@tmYLLabelsOn   = False    ; no y left tm labels
 ;res@tiYAxisOn      = False    ; no y-axis title
  res@gsnLeftString  = "IPR ( cm~S~-3~N~s~S~-1~N~ ) mean"       ;  -~S~222~N~Rn
  plot(2) = gsn_csm_contour_map_ce(wks,prnb,res) 

  ;; line over color 

  res@cnLinesOn       = False    ; turn off contour lines
  res@cnFillOn        = True     ; turn on color
  res@cnLineThicknessF= 2.       ; thicker lines 
  res@cnMinLevelValF  =    0.      ; set min contour level
  res@cnMaxLevelValF  =   16.      ; set max contour level
  res@cnLevelSpacingF =    2.
  res@lbLabelStride   =    2       ; every other label bar label
  res@gsnLeftString  = "IPR ( cm~S~-3~N~s~S~-1~N~ ) 90~S~th~N~ percentile"       ;  -~S~222~N~Rn
  plot(3) = gsn_csm_contour_map_ce(wks,prmb,res)           ; draw second plot
 ;overlay(plot(1),plotover)                              ; now over lay plots

 ;delete(plotover) 

;;----------------------------------------------------------------------
;; 3rd panel 
;;----------------------------------------------------------------------

  res@gsnRightString  = "JJA"
  res@cnLinesOn       = False    ; turn off contour lines
  res@cnFillOn        = True     ; turn on color
  res@cnMinLevelValF  =    0.    ; set min contour level
  res@cnMaxLevelValF  =   16.      ; set max contour level
  res@cnLevelSpacingF = 2.
  res@lbLabelStride   = 2        ; every other label bar label
 ;res@tmYLLabelsOn   = False     ; no y left tm labels
 ;res@tiYAxisOn      = False     ; no y-axis title
  res@gsnLeftString  = "IPR ( cm~S~-3~N~s~S~-1~N~ ) mean"       ;  -~S~222~N~Rn
  plot(4) = gsn_csm_contour_map_ce(wks,prnc,res) 

  ;; line over color 

  res@cnLinesOn       = False    ; turn off contour lines
  res@cnFillOn        = True     ; turn on color
  res@cnLineThicknessF= 2.       ; thicker lines 
  res@cnMinLevelValF  =    0.      ; set min contour level
  res@cnMaxLevelValF  =   16.      ; set max contour level
  res@cnLevelSpacingF =    2.
  res@lbLabelStride   =    2       ; every other label bar label
  res@gsnLeftString  = "IPR ( cm~S~-3~N~s~S~-1~N~ ) 90~S~th~N~ percentile"       ;  -~S~222~N~Rn
  plot(5) = gsn_csm_contour_map_ce(wks,prmc,res)          ; draw second plot
 ;overlay(plot(2),plotover)                         ; now over lay plots


;;----------------------------------------------------------------------
;; panel plot 
;;----------------------------------------------------------------------

  resp                    = True   ; panel mods desired
  resp@gsnFrame           = False  ; don't advance frame yet
 ;resp@gsnPanelBottom     = 0.0    ; space for label bar
 ;resp@gsnPanelTop        = 0.5    ; only panel on lower half of page
 ;resp@gsnPanelLabelBar   = True   ; common label bar
 ;resp@pmLabelBarWidthF   = 0.8    ; label bar width
 
  resp@gsnPanelYWhiteSpacePercent = 5
  resp@gsnPanelXWhiteSpacePercent = 5
 
  gsn_panel(wks,plot(0:5),(/3,2/),resp)

  frame(wks) 

  ;; output plot information 


end
 
          


